
# Analysis manuscript

rm(list = ls())

# This script contains the code to generate the figures and tables of the main manuscript.

#Set your path here

setwd("")

#Install and load required packages
library("pacman")
p_load(panelView, PanelMatch, fect, patchwork, texreg, estimatr, marginaleffects, grid, gridExtra, fixest, lubridate, forcats, 
       stringr, dplyr, purrr, readr, tidyr, tibble, tidyverse, ggspatial, ggplot2, sf) 

# Figure 1 ----------------------------------------------------------------

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")
shp <- st_read("replication/database/shp/mapa_municipios_ESPCAN.shp")
colnames(shp)[10] <- "mun_code"

map <- left_join(df, shp, by = "mun_code")
rm(df, shp)

map <- map |>  
  filter(year == "2011" | year == "2015" | year == "2016" | year == "2019")

map$depo_cat_te[map$depo_cat_te == "1_no_change"] <- "No change"
map$depo_cat_te[map$depo_cat_te == "2_decrease"] <- "Decrease"
map$depo_cat_te[map$depo_cat_te == "3_increase"] <- "Increase"

colours <-  c("#909497","#000000","#E5E7E9")
names(colours) <- c("No change", "Decrease", "Increase")

map |> 
  ggplot() + 
  geom_sf(aes(fill = depo_cat_te, geometry = geometry), size = NA, colour = "transparent", show.legend = T) + 
  coord_sf(datum = NA) +
  facet_wrap(~year, ncol = 2) +
  scale_fill_manual(values = colours, name = "Population change") +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 20),
        legend.position = "bottom",
        legend.title = element_text(size = 18),
        legend.text = element_text(size = 16))

ggsave(filename = "replication/outcomes/images/Figure1.png", width = 15, height = 10)



# Table 1 -----------------------------------------------------------------

rm(list = ls())
library("fixest")

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/long_time.RData")

m1 <- feols(c(psoe, pp, nswp ) ~
              sw(depo_cat_te) | mun_code + year,
            data = df,
            cluster = c("region", "mun_code" , "year")) 

setFixest_dict(c(psoe = "PSOE", pp = "PP", nswp="NSWP",
                 mun_code = "Municipality", 
                 year = "Year", 
                 region = "Region",
                 "depo_cat_te2_decrease"="Decrease",
                 "depo_cat_te3_increase"="Increase"))

etable(m1,style.df = style.df(depvar.title = "", fixef.title = "", 
                              fixef.suffix = " FE", yesNo = "Yes"), 
       tex = TRUE, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
       order = c("Depopulation"), view = T,  replace = T, file = "replication/outcomes/tables/Table1.tex")


# Table 2 -----------------------------------------------------------------


rm(list = ls())

loadRData <- function(fileName){
  #loads an RData file, and returns it
  load(fileName)
  get(ls()[ls() != "fileName"])
}

df <- loadRData("replication/database/short_time.RData")

setFixest_dict(c(psoe = "PSOE", pp = "PP", podemos = "UP", cs = "Cs", 
                 nswp= "PANES",  es ="ES", vox = "VOX", 
                 "depo_cat_te2_decrease"="Decrease",
                 "depo_cat_te3_increase"="Increase",
                 agriculture_workers = '% workers agriculture', 
                 services_workers = '% workers services',
                 population_log = "(Log) Population", 
                 mun_code = "Municipality", 
                 year = "Year",
                 region = "Region",
                 var_older_te = "$\\triangle{}$% 60 y/o or older",
                 var_younger_te = "$\\triangle{}$% 16 y/o or younger",
                 unemployed = "% unemployed"))

m1 <- feols(c(psoe, pp, podemos, cs, nswp, es, vox ) ~ 
              sw(depo_cat_te) +
              agriculture_workers + services_workers + unemployed + 
              var_older_te + var_younger_te + population_log | mun_code + year,  
            data = df,
            cluster = c("region", "mun_code" , "year")) 

etable(m1,style.df = style.df(depvar.title = "", fixef.title = "", 
                              fixef.suffix = " FE", yesNo = "Yes"), tex = TRUE, signif.code = c("***"=0.01, "**"=0.05, "*"=0.10),
       order = c("Depopulation"), view = T,  replace = T,file = "replication/outcomes/tables/Table2.tex")


# Figure 2 ----------------------------------------------------------------

modelslog <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelslog[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~depo_cat_te*population_log +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te | mun_code + year")),
          cluster = c("mun_code"))
  }

modelslog <- lapply(modelslog, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                    condition = c("population_log")) + 
                      geom_rug(data = df, aes(x = population_log),  alpha = 0.2) +
                      ylab("Marginal effects") + xlab("") +
                      ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, 
                                                                linetype="dashed", color = "blue") + theme_bw())

psoe <-  modelslog$psoe + ggtitle("PSOE")
pp <- modelslog$pp + ggtitle("PP")
podemos <- modelslog$podemos + ggtitle("Podemos")
cs <- modelslog$cs + ggtitle("Cs") 
nswp <- modelslog$nswp + ggtitle("NSWP")
es <- modelslog$es + ggtitle("ES")
vox <-  modelslog$vox + ggtitle("Vox")

bottom <- textGrob("(Log) Population", rot = 0, gp = gpar(fontsize = 20))
g_sizemun <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 3, nrow = 3, bottom = bottom)


ggsave(g_sizemun, filename = "replication/outcomes/images/Figure2.png", width = 15, height = 10)


# Figure 3 ----------------------------------------------------------------

modelsage <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelsage[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*var_older_te + 
                        agriculture_workers + services_workers + unemployed + population_log + 
                         var_younger_te | mun_code + year")),
          cluster = c( "mun_code" ))
}

modelsage <- lapply(modelsage, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                    condition = c("var_older_te")) + 
                      geom_rug(data = df, aes(x = var_older_te),  alpha = 0.2) +
                      ylab("Marginal effects") + xlab("") +
                      ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, 
                                                                linetype="dashed", color = "blue") + theme_bw())

psoe <-  modelsage$psoe + ggtitle("PSOE")
pp <- modelsage$pp + ggtitle("PP")
podemos <- modelsage$podemos + ggtitle("Podemos")
cs <- modelsage$cs + ggtitle("Cs") 
nswp <- modelsage$nswp + ggtitle("NSWP")
es <- modelsage$es + ggtitle("ES")
vox <-  modelsage$vox + ggtitle("Vox")

bottom <- textGrob("\u0394 % 60 y/o or older", rot = 0, gp = gpar(fontsize = 20))
g_old <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 3, nrow = 3, bottom = bottom)

ggsave(g_old, filename = "replication/outcomes/images/Figure3.png", width = 15, height = 10)


# Figure 4 ----------------------------------------------------------------


modelspublic <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  modelspublic[[i]] <- df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*public_services +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log  |  mun_code + year")),
          cluster = c( "mun_code" ))
}

modelspublic <- lapply(modelspublic, function(x) plot_cme(x, variables = "depo_cat_te", 
                                                          condition = c( "public_services")) + 
                         ylab("Marginal effects") + xlab("") +
                         ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, linetype="dashed", color = "blue") + theme_bw())

psoe <-  modelspublic$psoe + ggtitle("PSOE")
pp <- modelspublic$pp + ggtitle("PP")
podemos <- modelspublic$podemos + ggtitle("Podemos")
cs <- modelspublic$cs + ggtitle("Cs") 
nswp <- modelspublic$nswp + ggtitle("NSWP")
es <- modelspublic$es + ggtitle("ES")
vox <-  modelspublic$vox + ggtitle("Vox")

bottom <- textGrob("\u0394 Public services (education and health centres)", rot = 0, gp = gpar(fontsize = 20))
g_publicservices <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 3, nrow = 3, bottom = bottom)

ggsave(g_publicservices, filename = "replication/outcomes/images/Figure4.png", width = 15, height = 10)


# Figure 5 ----------------------------------------------------------------

models_services <- list()
for (i in c("psoe", "pp", "podemos", "cs", "nswp", "es", "vox")){
  models_services[[i]] <-  df |> 
    filter(depo_cat_te %in% c("1_no_change", "2_decrease") ) |> 
    feols(as.formula(paste0(i, " ~ depo_cat_te*amenities +     
                         agriculture_workers + services_workers + unemployed + var_older_te + 
                         var_younger_te + population_log  |  mun_code + year")),
          cluster = c( "mun_code" ))
}

models_services <- lapply(models_services, function(x) plot_cme(x, variables = "depo_cat_te",
                                                                condition = c("amenities")) + ylab("Marginal effects") 
                          + xlab("") +
                            ggtitle(x$call$formula[[2]]) + geom_hline(yintercept = 0, linetype="dashed", color = "blue") + theme_bw())

psoe <-  models_services$psoe + ggtitle("PSOE")
pp <- models_services$pp + ggtitle("PP")
podemos <- models_services$podemos + ggtitle("Podemos")
cs <- models_services$cs + ggtitle("Cs") 
nswp <- models_services$nswp + ggtitle("NSWP")
es <- models_services$es + ggtitle("ES")
vox <-  models_services$vox + ggtitle("Vox")

bottom <- textGrob("\u0394 in the presence of amenities indicator", rot = 0, gp = gpar(fontsize = 20))
g_urbanfabric <- grid.arrange(psoe, pp, podemos, cs, nswp, es, vox, ncol = 3, nrow = 3, bottom = bottom)

ggsave(g_urbanfabric, filename = "replication/outcomes/images/Figure5.png", width = 15, height = 10)


# Figure 6 ----------------------------------------------------------------

shp <- st_read("replication/database/shp/mapa_municipios_ESPCAN.shp")
colnames(shp)[10] <- "mun_code"

map <- left_join(df, shp, by = "mun_code")
rm(df, shp)

map <- map |>  
  filter(year == "2011" | year == "2015" | year == "2016" | year == "2019")

colour <- c("#909497", "#000000")
names(colour) <- c("No", "Yes")

map |> 
  ggplot() + 
  geom_sf(aes(fill = depo_risk, geometry = geometry), size = NA, colour = "transparent") + 
  coord_sf(datum = NA) +
  facet_wrap(~year, ncol = 2) +
  scale_fill_manual(values = colour, name = "Depopulation risk") +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold", size = 20),
        legend.position = "bottom",
        legend.title = element_text(size = 18),
        legend.text = element_text(size = 16))

ggsave(file = "replication/outcomes/images/Figure6.png", height = 10, width = 15)


# Figure 7 ----------------------------------------------------------------

df <- loadRData("replication/database/short_time.RData")
names(df)

df$risknum <- NA
df$risknum[df$depo_risk=="Yes"]  <- 1
df$risknum[df$depo_risk=="No"]  <- 0

df$coden <- as.integer(df$mun_code)
df$yearn <- as.integer(df$year)

out.psoe <- df  %>%
  fect(psoe ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"), 
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "PSOE (ATT)", ylab = "Effect of depopulation risk on support for the PSOE", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

out.pp <- df  %>%
  fect(pp ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"), 
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "PP (ATT)", ylab = "Effect of depopulation risk on support for the PP", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

out.podemos <- df  %>%
  fect(podemos ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"), 
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "Podemos (ATT)", ylab = "Effect of depopulation risk on support for Podemos", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

out.cs <- df  %>%
  fect(cs ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"), 
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "Cs (ATT)", ylab = "Effect of depopulation risk on support for Cs", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

out.nswp <- df  %>%
  fect(nswp ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"), 
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "NSWP (ATT)", ylab = "Effect of depopulation risk on support for NSWP", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

out.es <- df  %>%
  fect(es ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"),  
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "ES (ATT)", ylab = "Effect of depopulation risk on support for ES", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

out.vox <- df  %>%   
  fect(vox ~ risknum  +   
         agriculture_workers + services_workers + var_older_te + 
         var_younger_te + unemployed + population_log, data = ., index = c("coden","yearn"), 
       method = "fe", force = "two-way", na.rm = T, se = TRUE) %>%  
  plot(main = "Vox (ATT)", ylab = "Effect of depopulation risk on support for Vox", 
       cex.main = 0.8, cex.lab = 0.8, cex.axis = 0.8, stats = "F.p", ylim = c(-4,4))

g_did <- out.psoe + ylab("Effect of depopulation risk") + out.pp + ylab("Effect of depopulation risk") +
  out.podemos + ylab("Effect of depopulation risk") + out.cs + ylab("Effect of depopulation risk") +
  out.nswp + ylab("Effect of depopulation risk") + out.es + ylab("Effect of depopulation risk") +
  out.vox + ylab("Effect of depopulation risk")

ggsave(g_did, file = "replication/outcomes/images/Figure7.png", width = 15, height = 10)


